lpf = tf(30,[1 30]);

% Connect low-pass filter to first output of Px
Pxdes = append(lpf,1) * Px;
set(Pxdes,'outputn',{'x-gap*' 'x-force'})

% Design the state-feedback gain using LQRY and q=1, r=1e-4
kx = lqry(Pxdes(1,1),1,1e-4);

dt = 0.01;
t = 0:dt:50;   % time samples

% Generate unit-covariance driving noise wx = [w-ex;w-ix].
% Equivalent discrete covariance is 1/dt
wx = sqrt(1/dt) * randn(2,length(t));

lsim(Px(1,2:3),':',clx(1,2:3),'-',wx,t)
% Specify model components
Hy = tf(7.8e8,[1 71 88^2],'inputn','u-y') 
Fiy = tf(2e4,[1 0.05],'inputn','w-iy') 
Fey = tf([1e5 0],[1 0.19 9.4^2],'inputn','w-ey')
gy = 0.5e-6 % force-to-gap gain

% Build open-loop model
Py = append([ss(Hy) Fey],Fiy);
Py = [-gy gy;1 1] * Py;
set(Py,'outputn',{'y-gap' 'y-force'});

% State-feedback gain design
Pydes = append(lpf,1) * Py; % Add low-freq. weigthing
set(Pydes,'outputn',{'y-gap*' 'y-force'});
ky = lqry(Pydes(1,1),1,1e-4);

% Kalman estimator design
esty = kalman(Pydes(2,:),eye(2),1e3);

% Form SISO LQG regulator for y-axis and close the loop
Regy = lqgreg(esty,ky);
cly = feedback(Py,Regy,1,2,+1);
% Compare the open- and closed-loop response to the white noise input disturbances.

dt = 0.01;
t = 0:dt:50;
wy = sqrt(1/dt) * randn(2,length(t));

lsim(Py(1,2:3),':',cly(1,2:3),'-',wy,t);